function edge_mark = mark_edge_longest(V, T, E, TE, ind)
% function edge_mark = mark_mesh_edge(V, T, E, TE, ind)
% mark refine edges, 1 means to do refinement later
% we note that this is important for refine results, you should be careful!
% assume that ind is a vector with the same length of T
%  ind(i) == 1 means triangle i need to be refine
%  otherwise , ind(i) = 0;

% reserve memory for output value

% bdr = find(ET(:,2)==0);
% node_flag = zeros(nv,1);
% node_flag(E(bdr,1)) = 1;
% node_flag(E(bdr,2)) = 1;

% degree = full(sum(VT,2));

[len_max_idx, len_max] = longest_edge(V,T);

edge_mark = zeros(size(E,1),1);

ind = find(ind==1);  % make unique in the list

% rule one: only mark longest edge for bad quality triangles
%  1->4 refine can prevent being worth further!
% ratio = cal_ratio(V, T(ind,:));
% is_bad = ratio > 1.5;
% T_good = ind(~is_bad);      % 1->4 for good quality tris
% edge_mark(TE(T_good,:)) = 1;
%edge_mark(TE(ind,:)) = 1;

%% only mark the longest edge for each mesh
while ~isempty(ind)
    for kk = 1:length(ind)
        tri = ind(kk);
        edge_mark(TE(tri, len_max_idx(tri))) = 1;  % mark the longest side of triangle kk
    end
    ne = edge_mark(TE(:,1)) + edge_mark(TE(:,2)) + edge_mark(TE(:,3));
    % try to find our whether the longest edge is refined
    flag = zeros(size(T,1),1);
    for tri = 1:size(T,1)
        flag(tri) = (edge_mark(TE(tri,len_max_idx(tri)))==0); % longest have not been refined!
    end
    
    ind = find(ne > 0  & (flag==1));
        
end

% is_bad = 1:length(ind);  % regard all elements are bad temporally
% 
% T_bad = ind(is_bad);  % there are two categries: with one large angle  and  one small angle
% cos_angle = cal_angle_cos(V, T, T_bad);
% is_angle_small = cos_angle > cos(pi/6); 
% is_angle_large = cos_angle < cos(2*pi/3); 
% 
% % with one large angle, refine the opposite edge
% type1 = sum(is_angle_large,2) == 1;  T1 = T_bad(type1); 
% T11 = T1(is_angle_large(type1,1) == 1); edge_mark(TE(T11,1)) = 1;
% T21 = T1(is_angle_large(type1,2) == 1); edge_mark(TE(T21,2)) = 1;
% T31 = T1(is_angle_large(type1,3) == 1); edge_mark(TE(T31,3)) = 1;
% 
% % with only one small angle:
% type2 = sum(is_angle_small,2) == 1; T2 = T_bad(type2); 
% T11 = T2(is_angle_small(type2,1) == 1); edge_mark(TE(T11,[2 3])) = 1;
% T21 = T2(is_angle_small(type2,2) == 1); edge_mark(TE(T21,[3 1])) = 1;
% T31 = T2(is_angle_small(type2,3) == 1); edge_mark(TE(T31,[1 2])) = 1;
% 
% type0 = (~type1) | (~type2);  % do 1->4 refine for not-so-bad triangles
% edge_mark(TE(T_bad(type0),:)) = 1;  

% % Most simple require: write down the edges of marked triangles
% n_eg_marked = edge_mark(TE(:,1)) + edge_mark(TE(:,2)) + edge_mark(TE(:,3));
% ind = find(n_eg_marked == 2);
% while (~isempty(ind))
%     
%     edge_mark(TE(ind,:)) = 1;
%     
%     n_eg_marked = edge_mark(TE(:,1)) + edge_mark(TE(:,2)) + edge_mark(TE(:,3));
%     ind = find(n_eg_marked == 2);
% end
% 

% 3. if boundary element is marked, mark all edges
% post process the edges' mark according to the edge length
% the boundary edge in this case is not refined 
% refineTFlag2 = zeros(size(T,1),1);
% refineTFlag2(ET((ET(:,2) == 0),1)) = 1;
% bnd_tri_mark = (refineTFlag2 == 1) & (refineTFlag==1);
% edge_mark(TE(bnd_tri_mark,:)) = 1;


% 5. treate boundary
% stop = 0;
% while (~stop)
%     n_mark = edge_mark(TE(:,1)) + edge_mark(TE(:,2))  + edge_mark(TE(:,3)) ;
%     idx12 = find(n_mark > 0); 
%     t11 = idx12((degree(T(idx12,1)) >= 8) + (edge_mark(TE(idx12,1))==1) == 2);  
%     edge_mark(TE(t11,1)) = 0;
%     t12 = idx12((degree(T(idx12,2)) >= 8) + (edge_mark(TE(idx12,2))==1) == 2);  
%     edge_mark(TE(t12,2)) = 0;
%     t13 = idx12((degree(T(idx12,3)) >= 8) + (edge_mark(TE(idx12,3))==1) == 2);  
%     edge_mark(TE(t13,3)) = 0;
%     if(length(t11) + length(t12) + length(t13) == 0)
%         stop = 1;
%     end
% end

% % 10. if the edge at the boundary and its length is large
% x21 = V(E(idx,2),1) - V(E(idx,1),1);
% y21 = V(E(idx,2),2) - V(E(idx,1),2);
% len = x21.*x21 + y21.*y21;
% bnd_mark = idx(len > 3*mean(len));
% edge_mark(bnd_mark) = 1;   % remark them
end



function cos_angle = cal_angle_cos(V,T,ind)
%function cos_angle = cal_angle_cos(V,T,ind)
% calculate all the consine of inner angle of triangles listed in ind
%
x21 = V(T(ind,2),1) - V(T(ind,1),1);
y21 = V(T(ind,2),2) - V(T(ind,1),2);
norm21 = sqrt(x21.*x21 + y21.*y21);
x32 = V(T(ind,3),1) - V(T(ind,2),1);
y32 = V(T(ind,3),2) - V(T(ind,2),2);
norm32 = sqrt(x32.*x32 + y32.*y32);
x13 = V(T(ind,1),1) - V(T(ind,3),1);
y13 = V(T(ind,1),2) - V(T(ind,3),2);
norm13 = sqrt(x13.*x13 + y13.*y13);
%
cos_angle = [ (- x13.*x21 - y13.*y21)./(norm13.*norm21), ...
                     (- x21.*x32 - y21.*y32)./(norm21.*norm32), ...
                     (- x32.*x13 - y32.*y13)./(norm32.*norm13)];
%
end

function [len_max_idx, len_max] = longest_edge(V,T)
% function [len_max, len] = longest_edge(V,T)
% After calling, len(:,1) = square of the longest edges and len_max is the local index of longest
x12 = V(T(:,1),1) - V(T(:,2),1);
y12 = V(T(:,1),2) - V(T(:,2),2);
x13 = V(T(:,1),1) - V(T(:,3),1);
y13 = V(T(:,1),2) - V(T(:,3),2);
x23 = V(T(:,2),1) - V(T(:,3),1);
y23 = V(T(:,2),2) - V(T(:,3),2);
len = [(x23.*x23 + y23.*y23),  (x13.*x13 + y13.*y13), (x12.*x12 + y12.*y12)];

% the first edge is current longest
len_max = len(:,1);  
len_max_idx = ones(size(len,1),1);  

% the second one is longer
idx2 = (len(:,2) > len_max);  
len_max_idx(idx2) = 2;  
len_max(idx2) = len(idx2,2);  

% the second one is longer
idx3 = (len(:,3) > len_max);  
len_max_idx(idx3) = 3;  
len_max(idx3) = len(idx3,3);  
end
